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ABSTRACT 

The purpose of this investigation is to develop a new model for the acceleration 
of the fast solar wind by nonlinear, time-dependent multidimensional MHD 
simulations of waves in solar coronal holes. Computational studies indicate 
that nonlinear waves are generated in coronal holes by torsional Alfven waves. 
These waves in addition to thermal conduction may contribute considerably to 
the accelerate the solar wind. Specific goals of this proposal are to investigate 
the generation of nonlinear solitary-like waves and their effect on solar wind 
acceleration by numerical 2.5D MHD simulation of coronal holes with a broad 
range of plasma and wave parameters; to study the effect of random disturbances 
at the base of a solar coronal hole on the fast solar wind acceleration with a 
more advanced 2.5D MHD model and to compare the results with the available 
observations; to extend the study to a full 3D MHD simulation of fast solar wind 
acceleration with a more realistic model of a coronal hole and solar boundary 
conditions. The ultimate goal of the three year study is to model the fast solar 
wind in a coronal hole, based on realistic boundary conditions in a coronal hole 
near the Sun, and the coronal hole structure (i.e., density, temperature, and 
magnetic field geometry) that will become available from the recently launched 
SOHO spacecraft. We report on our study of the effects of slow magnetosonic 
waves recently detected by EIT on the solar wind. 


1. Introduction 

To account for the observed properties of the solar wind Alfven waves were suggested as 
a source of momentum and heat, and studied in the linear regime (e.g., Alazraki & Courtier 
1971, Belcher 1971, Hollweg 1973, Jacques 1977, Davila 1985, Davila 1987, Ofman h Davila 
1995). The WKB approximation is the standard approach to study the effect of the Alfven 
waves on the solar wind. Recent linear models include many realistic features and consider 
the propagation of three-dimensional Alfven in a stratified, thermally conductive solar wind 
with WKB and non WKB fluctuations. 

In order to consider fully the solar wind acceleration by transverse Alfven waves, 
nonlinear effects must be included. This is necessary in order to couple the transverse 
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motions and the gradients of the magnetic pressure into the longitudinal motions. Nonlinear 
two- fluid one-dimensional MHD solar wind model with WKB Alfven waves were developed 
and compared favorably with observations (e.g., Esser et al. 1997). However, with the 
WKB approach the Alfven waves are not treated self consistently with the solar wind. 
Nonlinear effects due to the Alfven waves in the solar wind, such as fast and slow waves 
were deduced from the second order analysis of the MHD equations (Lou 1993). Recently, 
self consistent nonlinear computations of the solar wind acceleration by Alfven waves were 
performed in one spatial dimension (Lau and Siregar 1996, Boynton k Torkelsson 1996) 
and in two dimensions (Ong et al. 1997, Ofman k Davila 1997a, Ofman k Davila 1998). 

Ofman and Davila (1997a, 1998) (hereafter, OD) investigated the self-consistent 
nonlinear effects due to Alfven waves in coronal holes via numerical solution of the nonlinear 
2.5-D (3-D with azimuthal symmetry) resistive MHD equations in a spherical geometry 
with an (r, #)-inhomogeneous atmosphere and the hydrostatic initial state. They found that 
solitary waves are driven by the Alfven waves in a stratified coronal hole. Ofman and Davila 
(199fb) suggested that the solitary waves may explain the apparent broadening of ion 
emission lines observed by the SOHO Ultraviolet Coronagraph Spectrometer (UVCS). The 
waves that OD found propagate in a highly inhomogeneous radially divergent, dispersive, 
and dissipative medium, and are more complex than the classical sound solitons. Therefore, 
we refer to them as large amplitude nonlinear (LAN) waves. 

Ofman, Nakariakov, k DeForest (1999) studied the propagation of slow magnetosonic 
waves in plumes using the inviscid ID linear wave equation, and the nonlinear 2D MHD 
code in spherical geometry and found good qualitative agreement with the observed spatial 
and temporal evolution of the fluctuations in plume brightness obtained by EIT. 

The compressive viscosity coefficient, as calculated using Braginskii’s theory (Braginskii 
1965) in the solar corona is many orders of magnitude larger than the shear viscosity 
coefficient (Porter, Klimchuk, k Sturock 1994, Ofman, Davila, k Steinolfson 1994) and 
should strongly damp the propagation of the slow magnetosonic waves. For typical coronal 
parameter, we find that the slow waves with frequencies in the few millihertz range will be 
dissipated within one solar radius above the photosphere. 

The contribution of the low frequency slow waves to heating depends on the amount 
of energy flux available in these waves. The prevailing view based on Doppler shift 
observations (Athay k White 1979) is that a slow wave energy flux of ~ 10 4 erg cm -2 s -1 
that reaches the corona is not sufficient to balance the coronal energy losses. However, 
Porter, Klimchuk, k Sturock claim that if the observed nonthermal Doppler broadening of 
20-60 km s 1 of coronal emission lines (see, for example recent measurements of nonthermal 
line widths with the SUMER instrument on SOHO reported by Warren et al. 1997, and 
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Chae, Schuhle. & Lemaire 1998) is attributed to slow waves, than the energy flux is greatly 
underestimated. The direct observation of slow magnetosonic waves in plumes by EIT (e.g., 
DeForest and Gurman 1998) lends some support to the latter view. 


2. Results 


2.1. Dissipative Slow Magnetosonic Waves in Plumes 


2.1.1. Dissipative MHD Model Equations 


The main goal of the present numerical simulation of the slow magnetosonic waves in 
plumes is to investigate the role of dissipation, plume structure, and nonlinearity on the 
propagation of the waves. For this purpose we model the polar plume in the r - 0 plane 
with a 2D MHD nonlinear code in spherical geometry, where r is along the solar radial 
direction (r = 1 is the coronal boundary), and 9 is the latitudinal direction. We include 
viscous (compressive) and resistive dissipation, and the polytropic energy equation. We set 
the polytropic index 7 equal to 1. We plan to use more realistic energy equation in our 
future studies. 


The normalized MHD equations are 

dp 


at + v ' {pv) = 0 ' 

| + (v .v)v = -^vp * 

ot p F r r 2 p p 

<9B 

_ = V x (v X B) + 5" l V 2 B, 

^ + v ■ VT = (7 - l)(S v +S r - TV ■ v), 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 


where n is the viscous stress tensor (Braginskii 1965), S v is the viscous heating term and S r 
is the resistive heating term. When 7 = 1 the source terms do not affect the temperature. 


The quantities in the MHD equations are normalized with r r/R Q , t t/r A , 
v — » v/v A , B — * B/Bo, p — > p ! poo ■ P — ► p/po- The physical parameters are the Lundquist 
number S = t t /t a , the resistive time scale r r = i-nR^/uc 2 (where u is the resistivity, 
and c is the speed of light), the Alfven time scale t a = Rq/v a , the Froude number 
F r = v\Rq/(GM q ) (where G is the gravitational constant, and M Q is the solar mass), 
the Euler number E u = p 0 / ( poov A ) = C 2 Jv\ (= (3/ 2), the sound speed C s - ^/ypo/poo, 
the Alfven speed v A = B 0 / y/4irp^. The values of the magnetic field B 0 , pressure p 0 , and 
density poo are taken at the solar boundary of the plumes. 
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Here, \vg are inteiested in the dissipation of longitudinal waves. Therefor, we consider 
the radial component of the viscosity 


d 2 v r f dv^ ^ vA f 2 d In 770 \ 

dr 2 \ dr r ) \r dr ) 

where the compressive viscosity coefficient is given by (Braginskii 1965) 


(v-n), = -|jjo 


(5) 


Vo = 0.96nkT p T p , ( 6 ) 

with the proton collision time given by 

P iy/nriXe* 

In the above, m p is the proton mass, T p is the proton temperature, e is the electron charge, 
and A ss 20 is the Coulomb logarithm. For typical conditions at the base of the corona with 
1-4 x 10 I\ and n = 5 x 10 8 cm 8 we get t p = 2.8 s. The collision time is two orders 
of magnitude shorter than the typical slow wave period observed by EIT. At 1.2/?© the 
density drops by a factor of ~ 4 (assuming isothermal static atmosphere). This will result 
in ion collision time of 11 s, about 55 times smaller than the wave period observed by EIT. 
Thus, the plasma is in the collisional regime in the region of interest, and Landau damping 
can be neglected there. The proton collision frequency becomes comparable to the wave 
period (600 s) at about 3 R®. 

Note that the compressive viscosity coefficient varies mainly with temperature and 
can be written as 770 = 7700 ^' 2 . In the present single-fluid MHD model we assume that 
that T = T p — T e , where T e is the electron temperature, and we use the normalization 
Vo — ► Vo/(RqVaPoo)- With this normalization we have 770 = fjC s /v A - 

The heating due to the compressive viscosity is given by 



We calculate and plot the radial heating rate 2-nr 2 S v for various cases discussed below. 


2.1.2. Linear Theory of Slow Magnetosonic Waves in Viscous Fluid 

Consider the linear ID dynamics of slow magnetosonic waves propagating along the 
vertical magnetic field in an isothermal atmosphere. The magnetic field is taken to be 
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radially divergent. The atmosphere is spherically stratified 
background density given by 


Po — Pooe 




by gravity with the equilibrium 


(9) 


where R e is the solar radius, and H is the scale height, defined as H = 2 k B TR%/(GM Q m n ) 
where k B is Boltzmann’s constant, T the temperature, G the gravitational constant, 

Mq the Sun’s mass and the hydrogen mass mu. The scale height can be estimated as 
//(10 8 cm) ~ 61T’(10 6 K) (the factor in front of T is 51 when Helium is taken into account). 


In spherical coordinates, the linear, spherical slow magnetosonic waves radially 
propagating strictly along the magnetic field, with perturbations of plasma density p and 
radial speed ty, are governed by the equations 


OP 

dt 


1 0 / 2 \ 

7 2 ~dr V PoVr ) = °’ 


(10) 


Po 


dv T 

dt 


+ C 2 ^ 

+ Cs dr 


4 d 

+ gP ~ 


1 S 2 

7^r (T ^ 


= 0 , 


(11) 


where g — GMq/t 2 is the gravity acceleration and rjo is the compressive viscosity coefficient, 
and C s is the constant sound speed (see, § 2.1.1 below). Perturbations of the magnetic field 
and transversal components of the plasma velocity are not excited by the linear radially 
propagating slow wave considered. 


Equations (10) and (11) can be combined into the wave equation 


d 2 v r 

Po ~dt^ 



1 & { 2 \ 

7Fr {r P0 ' V) 


Assuming solutions of the form f(r)e~ iult 
perturbation p gives 


7^r {r P ^ ~ 3 Po d~r 


'1 d 

( 2 0v r \ 

r 2 dr ' 

[ dt) J 


-u) 2 p 


Cl 0_ ( 2 dp\ 

r 2 dr \^ r dr ) 


dp .4 Cl d A p 
~ + '3 m ^d7 = 0 


(13) 


Since the viscosity is small we have used the ideal momentum equation to substitute for ^ 
in the viscous term. The viscosity is only important in regions of high gradients. Therefore, 
we have kept only the highest order derivative of p in the viscous term. The solution of the 
wave equation is of the form 


P ~ 


Rq 


e 


r 



(14) 


where k is complex, and z — r — R G . Expression (14) shows that the absolute amplitude of 
the density perturbations decreases with height, while the ratio of the density perturbations 
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amplitude (14) and the local equilibrium density value may grow provided that the viscosity 
is small ( Irn(k z ) < R e /2H near r = IRq), 


P o 



gik z—iuj t 


(15) 


In the isothermal inviscid atmosphere, the growth rate of the ratio p/p 0 depends on the 
temperature T only. In the presence of viscosity, the damping rate is subtracted from the 
growth rate (see the discussion at the end of § 2.1.3). 

The approximate local dispersion relation of the dissipative waves is obtained by 
substituting the above solution (14) into the wave equation (13) and keeping only the 
largest terms 


k ~ i± _ , .2770 a; 2 

~C S 8// 2 H + 3p 0 C*' 


(16) 


The real part of this dispersion relation agrees with the ideal dispersion relation for spherical 
sound waves (Torkelsson, & Boynton 1998). 


2.1.3. Nonlinear Slow Magnetosonic Waves 


To understand basic physical mechanism responsible for the evolution of slow 
magnetosonic waves in the static, isothermal, spherically stratified atmosphere, we 
investigate weakly nonlinear dynamics of the waves under the assumption 


p/po 1. 


(17) 


We restrict our attention to quadiatically nonlinear processes only. Adding quadratically 
nonlinear terms to the right hand sides of Eqs. (10) and (11), we rewrite the equations as 


where 


dp 1 d / 2 x 
dt + 7 dr v PoVt ) = Nu 


dV T 


>dp 


P °~dt + C » + 9P = N 2 , 


dr 


No = 


N ' = ( r >9 ■ 

9 , v 1 d ( 2 2 \ 4 a fi a / , x 

~dt {pVr) ~ 7*d~r ( r Po?V ) + 3 P °dr 7d~r ( r Vr ) 


(18) 

(19) 

( 20 ) 


( 21 ) 
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Fig- 1. — Dependence of the linear spherical slow wave amplitude on the distance from the 
Sun for four different normalized viscosities r). The dotted line corresponds to rj — 10~ 10 , the 
solid line to f) = 1CT 5 , the dashed line to rj = 1(T 4 and the dashed-dotted line to fj = 10 -3 . 
The initial wave period is 600 s. The atmosphere is isothermal with T = 1.4 x 10 6 K. The 
relative amplitude is given by p'(/?)[po(ft)p'(l)]~ 1 . 


It is convenient to put dissipative terms on the right hand side as well. 


Equations (18) and (19) can be combined into the weakly nonlinear wave equation 


^P_Cld_( 2 dp\_ dp 
dt 2 r 2 dr V dr) 9 dr 


dNi 

dt 


1 <9/o N 

( r ' V! ) • 


( 22 ) 


The right hand side terms N\ and iV 2 contain v r , which can be expressed through p 
using the linear equation (10). Equation (22) becomes Equation (13) in the linear limit. 


In the following, we consider the dynamics of the slow waves, assuming 

Vo 


X 

^<<1 


and 


C s p oqA 


« 1 , 


(23) 


where A is the wavelength, in addition to assumption (17). Conditions (23) correspond to 
the WKB approximation and the approximation of weak dissipation, respectively. These 
approximations allow us to consider slow modification of the waves by the weak nonlinearity 
and dissipation in the WKB approximation. The latter means that we neglect the wave 
reflection from the radial inhomogeneity. 
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Fig. 2. Dependence of the nonlinear spherical slow wave amplitude on the distance from 
the Sun for four different initial amplitudes A = p(l)/rh 0 {l). The solid line corresponds to 
A = 10~ 6 , the dotted line to A = 0.02, the dash-dotted line to ,4 = 0.08 and the dashed 
line to A = 0.16. The initial wave period is 600 s. The atmosphere is isothermal with 
T = 1.4 x 10 6 K. The normalized viscosity is 77 = 3.2 x 10“ 4 (i.e., “typical solar value”). The 
relative amplitude is given by /y(/?)[/Oo(i?)^(l)] -1 . 


Considering an outwardly propagating wave, we transform into the moving frame 


£ r C s t, Ft — er, (24) 

where the parameter e is small, representing the slow modification of the wave by the weak 
nonlinearity, stratification and dissipation. This allows us to apply the standard method 
of slowly varying amplitudes. Keeping terms of the first order of e only, we have from 
Eqs. (19), (20) and (21) we get 


V T 



N x = - 


Cj_d_ 2 

Po 


4r/ 0 C s d 2 p 

3p 0 d?'- 


and, finally, arriving at the evolutionary equation, 


(25) 


^P_ j. (I M 9(R) \ ^ 1 dp 2770 d 2 p n 

OR ^ 2C 2 ) P+ Po (R) P d$ 3C sPo (R) dp ~ °- 


(26) 


Equation (26) is a spherical analog of the Burgers equation for slow magnetosonic waves. 
When all three factors of the wave evolution (the stratification, nonlinearity and dissipation) 
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Fig. 3. Evolution of temporal dependences of initially harmonical nonlinear spherical 
Alfven waves with the distance from the Sun. The solid curve corresponds to the solution 
a f T? = IRq, the dotted at R = 1.15/?q, the dashed to R! = 1.3/?q and the dash-dotted to 
R' = 1.5 R q . The initial wave amplitude is .4 = 0.08. The initial wave period is 600 s. The 
atmosphere is isothermal with T = 1.4 x 10 6 K. The normalized viscosity is f) = 3.2 x 10 -4 . 
The relative amplitude is given by p'(R)/p 0 (R). 


modeled by the second, third and forth terms of Equation (26), respectively, are negligible, 
the wave form p is independent of R (p(R) = const). Consequently, the wave keeps its 
shape in the moving (with the speed C s ) frame of reference. When some of the evolution 
factors are taken into account, the wave form p depends on R , so the wave changes with the 
distance from the Sun. 


For the following consideration, it is convenient to introduce the dimensionless variables 


R' = R/Rq, P = pj Poo, H 1 = H/Rq, p'q = Po/poo, 

where H is the scale height. W ith variables (27), equation (26) can be re-written as 

<W_ , (J_ 1 \ 1 ,dp' fj d 2 p' 

dR’ \R' 2 H'R'i) P + po(R') P d? p 0 (R') ~ °’ 

where fj = 


(27) 


(28) 


When both nonlinearity and viscosity (the third and the forth terms) are negligible, 
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Equation (28) is easily integrated, giving us 


p'{R') = p'(l)M exp 


2H' 



(29) 


The linear solution (29) resembles solution (6) for a harmonic wave, but Eq. (29) is valid 
for waves of arbitrary initial shape as long as they satisfy conditions (23). 

Both the effects of stratification and viscosity can be taken into account for a harmonic 
wave 

p(R', £') = p'(R') exp (tfcf'), (30) 

where k is a dimensionless wave number. Substituting expression (30) into Equation (28) 
(where the third term, corresponding to nonlinearity, is neglected), we obtain 


P' ~ P'(!)4 7 exp 


x exp 


b (* - b 

-nk* J\x.v[k (i - -L) m 


Expressing the integral through the exponential integrals Ei(x), we have 


= p'(l)^ ex P 


1 ( 

1 Ml 

2 H' V 

Rf). 


X 


exp (f,e {l - «-.*(-*) + j-,e* [B (-T) - Ei (-jTj) 


( 31 ) 


(32) 


Expression (32) shows that there are two different scenarios of the wave propagation. 
When f)k is less than the critical value ( fjk 2 ) c , the relative amplitude of the wave, // j po 
initially grows, reaches a maximum, and then decays. For 77 k 2 > ( 7 ^ ) c , the wave decays 
from the very beginning. Dependence of the wave amplitude on distance from the Sun for 
different values of the viscosity are shown in Figure 1. The dotted curve corresponds to the 
dissipationless solution (29). When the medium has finite viscosity, the wave either reaches 
a maximum and then decay's or decays from R f — 1 Rq. To determine the critical value of 
the viscosity (rjk 2 ) c , we expand the function p f / po into a Taylor series near R! — 1: 


p'/po ~ l + (^7 - 1 - p^ 2 ) {R' - l). 

From the condition that the coefficient with R' - 1 changes sign, we obtain 


(Vk 2 ) c 


1 

2 IP 


(33) 


(34) 
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For a given viscosity, waves with wavelengths shorter than 


are decaying at R' = 1 R e 
viscosities higher than 


A < 27ti 


2H'fj 


, 1-2H' 

For a given wave period P, waves are decaying at R’ 


V > 


— Vf--d 

2itRq) V2 H' ) 


(35) 
1-R© for 


(36) 


Consequently, only sufficiently long wavelength slow r waves are propagating along 
plumes. However, if the relative wave amplitude grows with the distance from the 
Sun, nonlinear effects, described by the third term in Equation (28), come into play. 

The nonlinear generation of higher harmonics transfers the wave energy into the short 
wavelength part of the spectrum, where the strong dissipation takes place. Thus, waves 
of higher amplitude are subject to stronger dissipation. This is shown in Figures 2 and 
3. Figure 2 shows the growth of the wave amplitude for three different initial amplitudes. 
(The curves are normalized on the initial amplitudes). The waves of higher amplitudes 
reach their maxima quicker. In the vicinity of the maximum, the wave amplitude decays 
due to nonlinear dissipation. Nonlinear evolution of the wave shape is shown in Figure 3. 
Initially, the outward propagating wave grows and its shape steepens departing from a 
simple harmonic shape. Finally, as the wave propagates outward, its amplitude decreases 
and the waves returns to the initial harmonic shape. 

The analytical theory allows us to understand main features of short wavelength 
dynamics of slow magnetosonic waves in the radially stratified static atmosphere, taking 
into account effects of weak nonlinearity and dissipation. However, the analytical theory is 
based upon assumptions (17) and (23). In particular, the assumption that the wavelength 
is much less than the scale height, is not well justified for modeling of the slow waves in 
polar plumes. Also, effects of two-dimensional structuring as well as radial plasma flows, 
which affects the waves, are out of scopes of the analytical theory. Thus, in in the next 
section we present the results of more realistic quantitative modeling of the v r aves using the 
nonlinear dissipative ID, and 2D MHD model. 


2.1.4- 2D MHD Results 

The purpose of the 2D MHD simulations is to investigate the effects of the plume 
density and temperature structure on the propagation and the dissipation of the slow r 
magnetosonic v r aves. The 2D MHD simulations are initiated with the normalized number 
density and temperature profiles that represent a polar plume 
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Fig. 4. The ^-dependence of the initial normalized density and temperature of the plume. 


n o (0,r) — n m j n + (1 — n min )e n / 2 )M j no(r), 
T(9, r) = T max - (T max - T min ) 


(37) 

(38) 


where the angular plume half-width d = 0.04 Rad, and the normalized density of the 
surrounding coronal hole n m j n = 0.3, with the normalization rio(l) = 1 at the center of the 
plume. The normalized temperature in the center of the plume T min = 0.8 is lower than 
the normalized temperature of the surrounding corona T mSiX = 1. The 9 dependence of n 0 
and T is shown in Fig. 4. In the uniform magnetic field, such structure is a waveguide for 
trapped fast and slow magnetosonic waves. The initial radial dependence of v r and n 0 is 
given by Parker’s (1963) isothermal solar wind solution. 

The initial magnetic field B = BQ(r)e T is uniform in the ^-direction and varies as 

We use the following typical coronal hole parameters in our 2D run n 00 = 5 x 10 8 
cm -3 , B 0 = 10 G, v A = 975 km s~\ C 3 = 152 km s' 1 , r A = 718 s, and we set S = 10 6 , 
m = 5 x 10“ 5 , and T = 1.4 x 10 6 A". 


For the upper boundary of the plumes and at the ^-boundaries we use open boundary 
conditions. We allow for incoming characteristics at the solar boundary by using 
zero-order extrapolation for the variables (e.g., Steinolfson & Nakagawa 1976). The density 
perturbation n\(9,r = l,t) is driven harmonically with the frequency u = 7.5 and an 



Fig. 5. A snapshot at t — 8.55t,4 = 102.2 min of the spatial variations of the perturbed 
density n \ , v r , n, and S v due to the propagating slow magnetosonic wave in the model polar 
plume. 

amplitude of n d = 0.08: 

m{0, r = 1, t) = n d e~^ e ~ n/2)/< ^ 2 sinut. (39) 

The magnetosonic waves that are excited with the above driver are not coupled to Alfvenic 
fluctuations since d jd<p = 0. Thus v 0 and B ^ remain zero throughout the simulation in the 
plume. The slow wave is coupled to vg and Bg. However, the perturbations in vg and Bg 
remain small in our runs. 

We solve the 2D MHD equations using the fourth order Runge-Kutta method for 
temporal integration and fourth order spatial differencing with 600 grid points in the radial 
direction and 40 grid points across the half-width of the plume. We use symmetry boundary 
conditions at the center of the plume. In Figure 5 the spatial variations in the perturbed 
density m, v T , n, and S v at t = 8.55 t a = 102.2 min are shown. 

The spatial variation in n i and v T are nearly in phase, as expected from the linear 
theory of traveling waves in a dissipative medium. It is evident that the plume width affects 
the shape of the waves and leads to trapping of the waves. However, small amplitude 
magnetosonic waves can be detected outside the plume. This could be due to the energy 
contained in the tails of the Gaussian shape driver (Eq. (39)), small leakage of the wave 
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Fig. 6.— The temporal variations of m and v r in the plume at r = 1.20 and 9 = 1.57 due 
to the propagating magnetosonic wave. 

energy outside the plume, or evanescent tails of the trapped modes. The lower temperature 
inside the plume compared to the surrounding corona results in lower sound speed there. 
This leads to the bending of the wave fronts, and the corresponding bending of the viscous 
dissipation layers. 

The amplitude of the waves is small (m/n 0 < 1), and the solution is nearly linear 
near the solar surface. As the magnetosonic waves propagate outwards, the quantity m/n 0 
increases and the nonlinear steepening of the waves becomes apparent. The steepening 
of the waves is limited by the presence of the viscous dissipation, and no shocks can be 
foimed. However, the sharper gradients near the wave fronts lead to enhanced viscous 
dissipation there. The heating is distributed along the plume by the propagating waves up 
to a distance determined by the dissipation length of the weaves. 

In Figure 6 we show the temporal variation of v r , and m at r = 1.20i? o and 9 = 1.57 
(i.e., near the central axis of the plume). It is evident that the modeled density and velocity 
fluctuations of the traveling waves are in phase, as expected from the linear theory. The 
nonlinear steepening of the wave fronts and the resulting asymmetry in the shape of the 
waves is evident in the temporal evolution. The steepening is limited by the presence 
of compressive viscosity. Nonlinearity also modifies the initial hydrostatic background 
equilibrium and leads to changes in the wave amplitude and phase speed with time. 

The effects of nonradial plume divergence on the propagation and the dissipation of 
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Fig. 7. A snapshot at t — 5.30t 4 = 63.4 min of the spatial variations of the perturbed 
density fi \ , v r , n, and S v due to the propagating slow magnetosonic wave in the nonradial 
model polar plume with d = 0.01 at r = 1 8©. 

the slow magnetosonic waves are investigated by adding a quadrupole component to the 
radial background field: 


B = 



3 . 84(3 cos 2 0 — 1) 
2r 4 


e r 


384 sin 2d 
2H 


e<?, 


(40) 


wheie we set 80 0.1 and 84 — 0.6 to get and B(1,7t/2)=1 and an asymptotic divergence 

rate of 10. At r — 28© the cross sectional area of the plume diverged by a factor of about 
3. At r = 18© we set d = 0.01, and w r e model this divergence rate in the density structure 
of the plume by substituting d -> d(r) = d( 1) 8 r (r, 7 r/2) 1/2 r in Equations (37)-(38). B e , and 
the 6 dependence of 8 r is neglected in the divergence factor due to the small angular extent 
of the plume. The initial radial dependence of the density is taken to be locally hydrostatic 
with the temperature T(r,9). The initial state evolves to a self consistent solution during 
the 2D simulation. 


The spatial dependence of the solutions for the nonradial plume is shown in Figure 7 
with n 0 = 5 x 10 s , u> = 7.5, T{6 — 7t/ 2) = 1.12 x 10 6 K. The faster-than-radial expansion 
of the plume is evident in the spatial structure of the density. The main difference between 
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Fig. 8 — The temporal variations of m and v r in the plume at r = 1.20 and 6 = 1.57 (solid), 
6 = 1.58 (dotted), 9 = 1.59 (dashed), 9 = 1.60 (dash-dot), 9 = 1.61 (dash-dot-dot-dot) due 
to the propagating magnetosonic wave for the case in Figure 7. 

the radial and the nonradially divergent plume is in the morphology of the wave, S v , and 
the v r . The perturbations of the density are excited at the base of the plume and expand 
across the plume as they propagate upward. The density fluctuations are in phase with 
velocity fluctuations. The enhanced radial velocity of at the edges of the plume is due to 
the increase pressure gradient at the plume boundary. The viscous dissipation of the waves 
is more concentrated at the center of plume than in the radial case. 

The temporal evolution at r = 1.2 R Q at several locations across the plume and the 
coronal hole are shown in Figure 8. It is evident that the phase speed, and the amplitude 
of the waves \ aries considerably across the plume. At r = 1 . 27 ?q initial state evolves to 
a self-consistent solution within about 1.577,. The long term variations are due to the 
nonlinear evolution of the density and velocity in the plume and the surrounding coronal 
hole. 


2.2. Two Fluid Coronal Hole Model 

Neglecting electron inertia (m e < m p ; V = V p ; p = m p n), and relativistic effects 
(V < c), assuming quasi-neutrality (n = n e = n p ), neglecting minor ions, we get the 2-fluid 
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MHD equations 

% + V • (pV) = 0, (41) 

= -(Vb, + Vp p ) - 21hl + ij x B - V . n, (42) 

SB 1 1 Vn 

-— = -cV x E, E = —V x B + 77J + —J xB-^, (43) 

ot c nee ne 

V X B = — J, V e = V — — , (44) 

c ne v ; 

~of = ~7pPpV • V - V • Vp p + 5*, = -7 e p e V • Ve - Ve • Vp e + S e , (45) 

In the above equations we use the following normalization: r -> r/R s ; t -> t/r A \ V -> VjV A \ 
B _H " B/Bq\ p — * p/po; Pp -* Pp/Ppo! Pe — ► Pe/PeO- 5 = T r /r 4 the Lundquist number, 
r r = 47ra 2 /i/c 2 the resistive time scale, = a/V^ the Alfven time scale. In coronal holes 
for the relevant scales of 0(0.01 — 1/?q) the Hall term and the electron pressure terms are 
small in the Ohm’s law compared to the other terms, and V p « V e . The Ohmic heating 
term in the electron energy equation is S e = ( 7e - 1 )j 2 /S, Pi is the viscous stress tensor, 
and S v is the viscous ion heating term. 

To model the effect of random disturbances on the generation of waves we use a 
broad-band Alfven wave source at r = IR Q given by 


av 

dt 


+ (V- V)V 


0, r = 1) = -V d /V A , r F(t, 9 ) (46) 

100 

Q ) = Y 1 a isin{uit + r,-(0)) (47) 

i=i 

where a, no i 1 > no = 0. 115, u>, = u>\ -I- (i 1) Au, A u> (u>n — ui\ )/{N — 1), and r ( is a 
random phase. We use frequencies in the range of 1.4-6. 9 mHz. 


2.2.1. Results of the 2- fluid. Model 

In Figures 9-13 we show the results of the two-fluid run with the physical parameters 
7p = 11) 7e = 105, T p o = 2 x 10 6 K, T e o = 0.8 x 10 6 K, no = 10 8 cm -3 , So = 7 G, 

V A = 1527 km/s, Vd = 0.06. The Lundquist number is set to S = 10 4 . 

The spatial dependence of the velocity components and the density at t = 130r 4 = 16.5 
hrs are shown in Figure 9. The Alfvenic velocity component v # (bottom right panel) 
fluctuates out-of-phase with D 0 (see Figure 10) as the waves propagate in the coronal hole. 
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Fig. 9. The spatial dependence of the density and the velocities in the model coronal hole 
at t = 130t 4 = 16.5 hrs. 


The fluctuations in Vp were used to calculated the effective proton temperature discussed 
below. The v T (bottom left panel) reached about 400 km/s near r = 10 R Q and has not 
approached the asymptotic speed in this region. The fluctuations in v r in the central region 
of the coronal hole are due to the nonlinear compressional waves driven by Alfven waves. 
The vq remains small throughout the coronal hole. We have used the linearly polarized 
Alfven waves as the driver at the coronal boundary. Including other polarizations may 
produce larger Vg. 

In Figure 10 we show the corresponding magnetic field components and j 2 /p 
(proportional to the Ohmic heating per particle). It is evident that the largest gradients 
in Bp occur at the boundaries of the coronal hole, defined by the transition from low to 
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Fig. 10. The spatial dependence of the magnetic field components and j 2 / p in the model 
coronal hole at t = 130t 4 = 16.5 hrs. 

high density regions. The phase-mixing of the Alfvenic wave fronts leads to the generation 
of currents (see, also in this issue). The large currents in these regions lead to heating 
of the electrons via Ohmic dissipation, which is evident in the electron temperature (see 
Figure 11) below. The fluctuations in B d are small compared to B $ in accordance with the 
small fluctuations in v e . The variations in £ r /p 0 5 are due to the effects of the nonlinear 
compressional waves. 

The 2D distribution of the electron temperature in the model coronal hole is shown in 
the left panel in Figure 11. The electron temperature decreases from 0.8 x 10 6 K at the 
base of the corona to about 0.6 x 10 6 K at 10/? Q . Two regions of somewhat higher electron 
temperature are evident in the model coronal hole in the boundary between the low and 
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Fig. 11. — The spatial dependence of the electron and proton temperatures in the model 
coronal hole at t = 130^ = 16.5 hrs. 


the high density regions. The higher temperature is the result of the Ohmic heating due 
to phase-mixing of the Alfven waves in the inhomogeneous regions (see j 2 /p in Figure 10). 
The pioton tempeiatuie (right panel) is not affected by the Ohmic heating and decreases 
nearly monotonically from 2 x 10 6 K at r = 1 R Q to 0.9 x 10 6 at 10i? o . 

The radial dependence of the electron and proton temperatures and of v# and v r are 
shown in Figure 12. The fluctuations due to the nonlinear waves driven by the Alfven waves 
are evident in the v r , Tp , and T e . The nonlinear waves contribute to the acceleration of the 
solar wind in addition to the thermal gradient. 

In Figure 13 we show the temporal evolution of the velocity components (top panel) 
and the density perturbation (bottom panel) at r = 4 . 98 / 2 q at the center of the coronal 
hole. It is evident that the fluctuations due to the nonlinear waves in v r and p\ are in phase 
as would be expected for compressional longitudinal wave. The Alfvenic fluctuations of the 
v <t> are dominated by higher frequency than the nonlinear compressional wave. 

In Figure 14 we show the effect of v<j, on the line-width temperature. The dashed line 
show's the time-averaged proton temperature obtained form the model. This temperature 
is decreasing as expected from the Polytropic energy equation. The solid line shows the 
effective line-width temperature that incorporates the effect of waves on the simulated line 
width. The effective temperature calculated using T p , e// = T p + ^ >, is due to 

unresolved wave motions integrated over the line of sight in space and about 4 hours in 
time. It is evident that the effective line-width temperature increases with r in qualitative 
agreement with SOHO UVCS Ly-a observations Kohl et al. 1998. In the period covered by 
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Fig. 12 — The radial dependence of T p (solid line) and T e (dashed line) at the center of the 
model coronal hole are shown in the top panel. The v r (solid line) and the (dashed line) 
are shown in the bottom panel, 

this report we have further explored the parameter space of the solar wind with the 2.5-D 
two-fluid model. We have including the effects of self-consistent heating of the electrons 
by Ohmic dissipation, and heating of protons by viscous dissipation. We found that with 
millihertz range Alfven waves the solar wind will be heated beyond several solar radii, 
provide that strong phase mixing of the Alfven waves is present. 




3. Analysis 

Motivated by recent observations of slow magnetosonic waves in polar plumes 

by the EIT we investigated the propagation, and the dissipation of these waves. We 

derived the linear dispersion relation for these waves in the viscous plume plasma, and 
spherical geometry. We included the effects of weak quadratic nonlinearity and derived an 
evolutionary Burgers’ type equation for these waves with the main assumptions of weak 
nonlinearity, short wavelength (compared to the scale height), and no background solar 
wind flow. Numerical solution of the evolutionary equations show the dependence of the 
wave amplitude with height. To relax the above assumptions we used the ID and the 2D 
MHD equations to model the propagation, and the dissipation of the slow waves in plumes. 
The MHD models included compressive viscosity, self-consistent nonlinear effects, and 
background solar wind flow. The results of the simulations describe the dynamics and the 
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Temporal plot al r= 4.98164; theta= 1.57406 



Temporal plot at r= 4.98164; theta= 1.57406 



Fig. 13.— The temporal evolution of v r (solid line, top panel), (dashed line, top panel) 
and the perturbed density (bottom panel) in the center of the coronal hole. 


dissipation of the waves for a range of coronal parameters. 

The required energy flux to heat the corona (i.e., to balance the radiative and 
conductive losses) is about 10 5 ergs cm" 2 s -1 . The energy flux carried by the slow 
magnetosonic waves is given by p { ^-C s in the WKB approximation. For n = 5 x 10 8 
cm and T 1.4 x 10 K, the required wave amplitude to obtain sufficient energy flux 
to balance the losses is 40 km s h This wave amplitude is consistent with the observed 
nonthermal broadening of emission lines (e.g, Chae, Schuhle, & Lemaire 1998). The direct 
observation of slow magnetosonic waves in coronal plumes by DeForest and Gurman (1998) 
using the EIT instrument shows that these waves are present in coronal hole plumes. 

Slow magnetosonic waves with an amplitude of 40 km s -1 (about 27% of the sound 
speed) carry an order of magnitude more energy than the ~15 km s~ ] that were seen in EIT 
observations. The 40 km s _1 waves are in the nonlinear regime with significantly shorter 
dissipation length than the length predicted by the linear theory. We found numerically 
that due to nonlinear steepening of the wave fronts the dissipation length near the Sun is 
about 0.087?© for these waves, compared to 0T4.Ro for the 15 km s' 1 waves with 300 s 
period and solar value of compressive viscosity. 

The short dissipation length of high amplitude five-minute slow magnetosonic waves 
will make them hard to detect by the EIT instrument, since the signal would be dominated 
by the longer period, smaller amplitude waves that have longer dissipation length. In 
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Fig. 14.- Effective (i.e., including the effect of waves) and model line-of-sight time-averaged 
proton temperatures. 

addition, the available cadence of EIT does not allow observations of waves with periods 
below ~6 minutes. The UVCS WLC observations at 1.9/2© reported by Ofman et al. (1997, 
1998) are at higher cadence of ~ 90 s. Our simulations show that most of the waves in the 
millihertz frequency range would be dissipated by 1.9/2© for solar value of the compressive 
viscosity with T = 1.4 x 10 6 K. The detection of these waves at 1.9/2© possibly indicate that 
the effects of compressive viscosity are somewhat lower allowing the slow wave to propagate 
to 1.9/2©. This could be due to somewhat lower temperature, or higher density than the 
values used as “typical” solar values in our study. It is also possible that the slow waves are 
generated locally above 1/2©. 

We extended OD model by using more realistic parameters and source terms in the 
MHD equations. We improved the numerical technique and geometry of OD, replacing 
the hydrostatic initial state with Parker’s (1963) isothermal solar wind solution and with 
more complex initial magnetic field, and source terms. In particular, we found that the 
low frequency Alfven waves with a period of about an hour drive LAN waves in the radial 
flow velocity, which contributes significantly to the acceleration of the solar wind in coronal 
holes such that the fluctuating solar wind speed at 32/2© agrees with the observed fast solar 
wind speed of 500 - 900 km s" 1 (Grail et al. 1996). We found that in order to heat the 
solar wind to 10 K we need to deposit the heating close to the Sun by including strong 
phase mixing in solar wind inhomogeneities. We found that the dynamical evolution and 
the spatial variability of the LAN wave driven wind are in better qualitative agreement 


24 - 


with observations than the slowly varying, thermally driven and WKB Alfven wave solar 
wind models (for a review of recent models and observations see Winterhalter et al. 1996). 

We included a broad-band spectrum of the Alfven waves and investigated the effect of 
two-temperature electron-ion plasma of the solar wind on the signatures of the waves that 
would be expected from the SOHO UVCS observations. For this purpose we have solved 
the two-fluid MHD equations and found that the effective proton temperature in the model 
increases with height. The extension of our model into the two-fluid description with a 
broad band Alfven wave source and more realistic heating terms enables us to include more 
realistic electron and proton temperatures. 

The SOHO/UVCS observations suggest that T v > T e in coronal holes. We use the 
two-fluid 2.5D MHD model to investigate the self-consistent acceleration and heating of the 
two-temperature electron-proton plasma by a broad band Alfven wave source. We calculate 
the effective line-width proton temperature using the results of the 2.5D two-fluid model, 
an increase of the line-width proton temperature with heliocentric distance. This effect is in 
qualitative agreement with UVCS observations of the increase of the Ly-ct line width with 
height. 

Our coronal hole model suggest that the observed increase of the line-width proton 
temperature is due to temporally and spatially unresolved low-frequency wave motions of 
the protons. Thus, the MHD waves may contain a significant fraction of the solar wind 
energ\ in the lower corona, and therefore, may play an important role in the acceleration 
and heating of the fast solar wind. 


4. Recommended further action 

We plan to improve the energy equation in our model by relaxing the polytropic 
assumption and to investigate the effects of a realistic heat conduction term. Due to the 
small time step required to correctly describe the thermal conduction with an explicit 
numerical method we are able to include the heat conduction terms in the ID version of our 
2D code. In order to eliminate this limitation and use the full 2D capabilities of our code 
we plan to use an implicit method for the energy equation. This will allow to include the 
heating terms due to waves more realistically in our model. We will also plan to extend the 
plume study to the multi-fluid model with the more realistic energy equations for proton 
and electrons, and with self consistent heating terms. We plan to further improve the 
execution speed of our code by continue the optimization for parallel execution on the Cray 
T3E. 
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In the third year of this project we plan to expand our model in two ways: (1) we plan 
to include a second ion fluid in the equations, and thus to develop a three-fluid model. 
This is needed in order to understand recent SOHO/UVCS observations of the solar corona 
in Oxygen emission at 1032A and in other minor ion emission; (2) we plan to expand our 
model into a full 3D description, which will allow to use more realistic model of a coronal 
hole and solar boundary conditions. For this purpose we will include the djdcfi terms in our 
single fluid model, and investigate various methods of solutions, and the effect of boundary 
conditions. Initially, we will develop the code on our SGI workstations. Next, we will 
investigate the timing and the efficiency of the solution on the Cray J90 and the Cray T3E. 
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